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Human control behavior is rarely completely stationary over time due to fatigue or 
loss of attention. In addition, there are many control tasks for which human operators 
need to adapt their control strategy to vehicle dynamics that vary in time. In previous 
studies on the identification of time-varying pilot control behavior wavelets were used 
to estimate the time-varying frequency response functions. However, the estimation of 
time-varying pilot model parameters was not considered. Estimating these parameters 
can be a valuable tool for the quantification of different aspects of human time-varying 
manual control. This paper presents two methods for the estimation of time-varying pilot 
model parameters, a two-step method using wavelets and a windowed maximum likelihood 
estimation method. The methods are evaluated using simulations of a closed-loop control 
task with time-varying pilot equalization and vehicle dynamics. Simulations are performed 
with and without remnant. Both methods give accurate results when no pilot remnant 
is present. The wavelet transform is very sensitive to measurement noise, resulting in 
inaccurate parameter estimates when considerable pilot remnant is present. Maximum 
likelihood estimation is less sensitive to pilot remnant, but cannot detect fast changes in 
pilot control behavior. 
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sinusoid amplitude, deg 

tracking error signal, deg 

probability density function, - 

target forcing function, deg 

Heaviside step function 

vehicle dynamics 

open-loop dynamics 

pilot visual response 

cost function value, - 

imaginary unit 

vehicle dynamics gain, - 

pilot remnant gain, - 

pilot visual gain, - 

sinusoid index 

Fisher information matrix 

number of points 

pilot remnant signal, deg 

forcing function frequency integer factor 

remnant/control signal power ratio, - 

Laplace variable 


T c 

vehicle dynamics time constant, s 

Ti 

pilot lead time constant, s 

T m 

measurement time, s 

T n 

pilot remnant time constant, s 

t 

time, s 

u 

pilot control signal, deg 

y 

controlled system output signal, deg 

w 

wavelet transform, - 

w 

Gaussian white noise signal, - 

X 

time sequence, - 

Symbols 

a 

line search parameter, - 

e 

prediction error, deg 

fit 

sinusoid phase shift, rad 

Fm 

phase margin, deg 

V 

nondimensional time, - 

A 

wavelet scale, s 

0 

parameter vector 
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a 

standard deviation 

w t 

sinusoid frequency, rad s _1 

T v 

pilot visual delay, s 

w 0 

Morlet wavelet parameter, - 

UJ 

frequency, rad s ~ 4 


wavelet function, - 

UJ C 

crossover frequency, rad s _1 

i>o 

mother wavelet function, - 

M nm 

pilot neuromuscular frequency, rad s ^ 1 

C nm 

pilot neuromuscular damping, 


I. Introduction 

Techniques for the identification and prediction of single- and multi-loop time-invariant human control 
behavior are well established . 1-4 In order for these techniques to provide accurate identification results, 
specific requirements are posed on the design of experiments in which human control behavior is estimated. 
These requirements ensure that the human controller behaves like a stationary control element in the control 
loop. However, in real-life control tasks, human control behavior is rarely completely stationary over time 
due to fatigue or loss of attention. In addition, there are many control tasks for which human operators 
need to adapt their control strategy to vehicle dynamics that vary in time. Examples are the manual control 
during mode transitions in tilt-rotor aircraft, and the manual control of rotorcraft with variable speed rotors 
or aircraft with reconfigurable flight control systems . 5 

In the past decades, significant research has been devoted to the modeling of time-variant pilot control 
behavior . 6 However, most of this research was focused on the development of complex mathematical models 
to describe time-varying behavior and not on the estimation of time-varying frequency responses or time- 
varying model parameters. An accurate estimate of these parameters allows for a quantification of the change 
in pilot control behavior over time and gives valuable insight into how a pilot adapts to time- varying controlled 
dynamics or environmental variables. Wavelet transforms have been used in many different applications - 
such as, geophysical applications - to analyze localized variations of power within a time series , 7,8 but have 
only recently been introduced into the field of systems and control. 

In recent years, work on the identification of time-varying pilot control behavior focused on the use of 
wavelet transforms to identify time-varying frequency response functions and derived parameters such as 
crossover frequencies and phase margins . 9,10 However, several difficulties arise when using wavelets. First, 
there are many types of wavelets to consider, each with their advantages and disadvantages in terms of 
energy localization, rise time, and smoothness in the time and frequency domains. In addition, each wavelet 
has one or more parameters that define its properties in the time-frequency plane. The type of wavelet and 
its parameters need to be selected by the researcher, making this technique not straightforward to use . 11 
Second, there are insufficient techniques to fully asses the accuracy (bias and variance) of the estimated pilot 
frequency responses. 

In previous studies on the identification of time-varying pilot control behavior using wavelets, the es- 
timation of time-varying pilot model parameters was not considered . 9, 10 When attempting to estimate 
time-varying model parameters using wavelets a second step is required, fitting the pilot model frequency 
response to the wavelet transform frequency response. This is similar to the two-step methods used for the 
estimation of time-invariant pilot model parameters . 2 3 A technique exists to estimate the model parameters 
from time-domain data in a single step using maximum likelihood estimation (MLE) , 4 significantly reducing 
bias and variance compared to the two-step methods. However, this technique has not been used for the 
estimation of time-varying parameters in the field of manual control. 

This paper presents the results of a study on the estimation of time-varying pilot model parameters to 
quantify time-varying human control behavior. Two methods are considered. The first is a two-step method 
using wavelets for the identification of the time-varying pilot frequency responses in the first step. The 
second method is the use of MLE to estimate the time-varying parameters directly from time-domain data 
in a single step. A simulation with time-varying pilot and system dynamics was set up to evaluate both 
methods. The paper is structured as follows. First, the pilot control behavior identification problem will be 
discussed in Sec. II. The parameter estimation methods using wavelets and MLE will be discussed in Sec. Ill 
and Sec. IV, respectively. Next, the simulation setup will be given in Sec. V, followed by the identification 
results in Sec. VI. Finally, a discussion and conclusions will be given in Sections VII and VIII, respectively. 
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II. Identification of Pilot Control Behavior 

The human operator is a non-linear biological system. However, in a continuous control task when 
trained properly and given constant environmental conditions - the operator’s manual control behavior 
can be described by a quasi-linear time-invariant model and a remnant signal that accounts for non-linear 
behavior. 1 This cybernetic approach to characterize human control behavior has been a powerful research 
tool to investigate the effects of different perceptual cues on pilot control behavior, 12 assess aircraft handling 
qualities, 13 and evaluate different control system designs. 14 

A single- loop compensatory control task is presented in Figure 1. In this control task, a pilot is actively 
controlling the vehicle dynamics H c , while following a target signal ft- The error signal e that is the difference 
between the target forcing function f t and the vehicle dynamics output y is presented on a compensatory 
display. This task setup allows for the identification of a single frequency response function H p that defines 
pilot control behavior. The output of this linear response function combined with a remnant signal n that 
accounts for non-linear behavior, defines the pilot’s control output u. 

To excite the combined pilot-vehicle system for the identification of pilot control behavior, it is common 
to use a forcing function f t that is defined as a multi-sine signal. These signals have power at distinct input 
frequencies that cover the frequency range of human control. To reduce the difficulty of the control task and 
to lower the probability of crossover regression, the number of input frequencies is usually limited to 12 and 
the power at higher input frequencies is reduced. 1 

Current methods, used to characterize pilot control behavior, are based on the concept that the human 
controller is a linear time-invariant control element in the control loop. However, human control behavior is 
rarely time invariant, due to factors like fatigue or loss of attention. This means the linear pilot frequency 
response function H p is a function of time. Furthermore, there are many control tasks where human operators 
need to adapt their control strategy to vehicle dynamics H c that change in time, for example, manual control 
during mode transitions in tilt-rotor aircraft or manual control during a failure of aircraft flight control 
system components. Finally, the visual properties of cockpit instruments or displays might change in time, 
for example, due to a change in environmental lighting conditions. 

When using linear time-invariant system identification techniques, small variations of pilot control be- 
havior in time are captured in the remnant signal. However, when these variations are more significant 
and are the subject of study, identification techniques that are capable of capturing time-variant system 
dynamics need to be used. In previous research on the identification of time-varying pilot control behav- 
ior, the wavelet transform was used to identify time-varying pilot frequency response functions and derived 
parameters, such as crossover frequencies and phase margins. 9, 10 However, the estimation of time-varying 
pilot model parameters was not considered. 

Two different techniques to estimate the time- varying pilot model parameters are evaluated in this paper. 
An overview of these techniques is given in Figure 2. The estimation of the time-varying model parameters 
using wavelets requires two steps. In the fist step, a time-varying frequency response is identified using the 
wavelet transform. In the second step, the model parameters are estimated at each time step by fitting a pilot 
model to the frequency response. Using MLE to estimate the time-varying pilot model parameters requires 
only one step. The parameters are estimated directly using the time-domain data. A sliding time window 
will be used to estimate the parameters at each time step. These two parameter estimation techniques are 
discussed further in the following sections. 
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Figure 2: Identification methods for the estimation of time-varying pilot model parameters from measured 
time traces. 


III. Wavelet Transform Parameter Estimation 

This section describes the estimation of time- varying model parameters using the wavelet transform. The 
wavelet transform is used to identify a time- varying frequency response in the first step, followed by a pilot 
model fit to estimate the model parameters in the second step. Both steps will be discussed in this section. 

A. Step 1: The Wavelet Transform 

The wavelet transform can be used to analyze a time series that contains time-variant power at different 
frequencies. Numerous good references exist that give a complete overview of the continuous and discrete 
wavelet transforms, the different types of wavelet functions, and a large number of examples in different 
applications. 7, 8 ’ 11 This section provides a brief overview of the continuous wavelet transform and the Morlet 
wavelet as used in this study. Notations and equations are adopted from Ref. 8. 

The continuous wavelet transform of a discrete time signal x n , with n = 0, 1, . . . , N — 1 and N the number 
of points in the time series, is defined as the convolution of x n with a wavelet function ip that is a scaled 
and translated version of a mother wavelet ipo: 


N - 1 

w( a ) = x ™r 

m — 0 


(m — n)St 
A 


(1) 


with St the sampling interval and A the wavelet scale. * indicates the complex conjugate. By varying 
the wavelet scale A and translating along the localized time index n, the magnitude information can be 
constructed as a function of scale and time. If complex wavelets are used, magnitude and phase information 
can be constructed. When the scale A increases, the wavelet becomes more spread out in time and takes 
only low-frequency features of the time sequence into account, and vice versa. 

The convolution in Eq. (1) is performed N times for each scale, which is a computational intensive 
procedure. The continuous wavelet transform can be approximated by performing the calculations in Fourier 
space using the Fourier transform, making the procedure considerably faster. Using the convolution theorem, 
the continuous wavelet transform can be approximated by: 


N —1 

W( A) « M* ( Aw fc) eiulknSt (2) 

fc= o 

where X and tk are the Fourier transforms of x and ip, respectively, k is the frequency index and Wfc the 
angular frequency defined by: 


Wfc 


2irk 

N6t 


( k < N/2) and w*, = 


2nk 


(. k > N/2) 


( 3 ) 
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(a) time domain 



t/ A, - 


(b) frequency domain 



\uj/(2ir), - 


Figure 3: The Morlet wavelet function (wo = 6 and A = 105t). 


Using Eq. (2), the continuous wavelet transform for a given scale can be calculated simultaneously for 
the entire time sequence with N points. The wavelet function *F in Eq. (2) needs to be normalized to have 
unit energy at each scale A to ensure that each wavelet transform is weighted only by the amplitude of X 
and not by the wavelet function itself. This allows for the wavelet transforms at each scale to be directly 
comparable. The normalization is defined by: 


/27ts\ 1/2 

T(Acc fc ) = ( T 0 (Aw fe ) 


(4) 


There are many types of functions that can be used as mother wavelets. Commonly used functions for 
the continuous wavelet transform are the Morlet, Paul, and Mexican hat. A valid wavelet function must have 
zero mean and unit energy |*Fo(w)| 2 du> = 1). Furthermore, the function must have compact support 

or sufficiently fast decay to obtain localization in both the time and frequency domains. 

One of the most widely used continuous wavelets is the complex Morlet wavelet, which consists of a wave 
signal modified by a Gaussian envelope. The Morlet wavelet is used in the remainder of this report, as it 
has the property that it minimizes the time-frequency localization. 10 Furthermore, because of the complex 
nature of this wavelet, it is able to detect both time-dependent amplitude and phase information at different 
frequencies, which is a requirement for system identification applications such as the identification of pilot 
control behavior. The Morlet wavelet is defined by: 


MV) = 7T -1//4 e iaJor? e - ’ 72 / 2 (5) 

with wo = 12 the nondimensional frequency, and ?/ = t/X a nondimensional time parameter. The real and 
imaginary parts of this wavelet function are depicted in Figure 3a. The Fourier transform of the Morlet 
wavelet function is given by Eq. (6), in which H{ui) is the Heaviside step function. The frequency response 
function of the Morlet wavelet function is depicted in Figure 3b. 

T 0 (sw) = 7r- 1 / 4 Ff(w)e- (A " J -“ o)2 / 2 (6) 

Using the wavelet transform defined in Eq. (2), it is necessary to choose a set of scales A. For a given type 
of wavelet function, a relationship between scales and frequencies can be derived analytically by substituting 
a cosine wave of a known frequency into Eq. (2) and deriving the scale at which the wavelet power spectrum 
reaches its maximum. Using this approach for the Morlet wavelet function, scales can be related to frequencies 
using the following expression: 


^ _ W 0 + y/2 + wg 
2w 

with w the frequency in radians per second. 

As the power of the forcing function is limited to certain input frequencies, the linear pilot response 
function can only be identified at these frequencies. This means the continuous wavelet transform only needs 
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to be calculated at the input frequencies of the forcing function. This can be accomplished by calculating 
the scales corresponding to the forcing function input frequencies using Eq. (7). Using these scales, the 
time-varying pilot frequency response function can be approximated by: 




W u ( A) 
W e ( A) 


( 8 ) 


where W u and W e are the wavelet transforms of u and e, respectively, defined by Eq. (2). For the Morlet 
wavelet function, the scales are related to the input frequencies of the forcing function by Eq. (7). The same 
approach can be taken to estimate the time-varying frequency response function of the controlled dynamics 
or the time-varying open-loop response H 0 i(jui,t). 


B. Step 2: Pilot Model Fit 

The parameters of a pilot model are estimated by fitting the model to the estimated time-varying frequency 
response function from Eq. (8) by minimizing a weighted least-squares criterion using a constrained gradient- 
based optimization. The cost function used in the current study is defined by: 


J(0) 


iv 

N t Z - 


H p (juj) - H p (juj, 0) 


k= 1 




(9) 


where N t is the number of input frequencies where the estimated pilot frequency response function is defined 
and 0 is the vector of parameters to be estimated. This criterion quantifies the difference between the 
estimated pilot frequency response function H p and the frequency response function of the pilot model H p . 
By repeating this procedure at every time step where the pilot frequency response function is defined, the 
change of the pilot model parameters in time can be determined. To reduce the computational effort needed 
to calculate the time- varying parameters, the number of calculations can be reduced by resampling the time- 
varying frequency response to a lower sampling rate in time. In the current study, the initial parameters for 
the optimization problem are the (to be estimated) simulated pilot model parameters. 


IV. Maximum Likelihood Parameter Estimation 


Maximum likelihood estimation has been used in previous research to estimate the parameters of time- 
invariant pilot models. 15 An extensive description of MLE for the estimation of pilot model parameters is 
given in Reference 4. This section gives a brief description of the MLE procedure adopted in this paper. 

MLE attempts to find an estimate 0 of the parameter vector 0 that maximizes the likelihood function. 
The likelihood function L (0) is defined as the joint conditional probability density function of the prediction 
error for m measurements of the pilot control signal u(k): 


L (0) = /(e(l), e(2), . . . , e(k ), . . . , e(m)|0) (10) 

The prediction error, indicated as e(k) in Eq. (10), is defined as the difference between the measured 
pilot control signal u(k) and the modeled pilot control signal u(k) at discrete instants. When the remnant 
is assumed to be an additive zero-mean Gaussian white noise signal, the conditional probability density 
function for one measurement of e(fc) is given by: 


f«k) |0) 


1 * 2 (fc) 

g 2 <t^ 

V 27ra n 


( 11 ) 


The set of parameters that maximizes the likelihood function is the maximum likelihood estimate of the 
parameter vector 0. For the MLE method it is common practice to minimize the negative natural logarithm 
of the likelihood function instead of maximizing L (0), as this results in a more straightforward optimization 
problem. When a global minimum of the negative log-likelihood is attained, the resulting parameter vector is 
the maximum likelihood estimate, indicated with 0 ml- For a single-output pilot model, combining Eq. (10) 
and Eq. (11) results into the following expression for the maximum likelihood parameter estimate: 
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0ml = argmin —In L (0) = argmin 

e e 


m 2 
~2~ n 


m 


k = 1 


( 12 ) 


The minimum of this function is found using an unconstrained gradient-based Gauss-Newton optimiza- 
tion. The iterative parameter update equation for the Gauss-Newton optimization is given by: 


0 (i + 1) = 0 (i) - a (i) Mqq(Q (i)) dL( g e (i)) (13) 

where a is the line-search parameter. This parameter, which typically varies between 0 and 1, is determined 
at each iteration before the actual parameter update to ensure the most rapid minimization of the likelihood 
function. The gradients of the likelihood function with respect to all parameters, dL/( 90, can be evaluated 
using the Jacobians of the pilot model state-space matrices with respect to the parameter vector 0. The 
Fisher information matrix, indicated with the symbol Mqq in Eq. (13), is given by: 


i” .fa,(k)\ 

e9 “^Sv^0 -) 


(14) 


The MLE optimization requires the pilot models to be written in output-error state-space form. These 
state-space representations are easily obtained by converting the transfer function models using the con- 
troller canonical form. For the conversion to a state-space representation, the pilot model time delays are 
approximated using fifth order Pade approximations. Due to these Pade approximations, the matrices of the 
state-space model contain coefficients that are highly nonlinear functions of the pilot model parameters. 

To estimate time- varying pilot model parameters in the current study, 
the MLE optimization is performed at every time step tj using a sliding 
time window of length At. This is visualized in Figure 4. Choosing a 
At that is too small will decrease the accuracy of estimated parameters 
related to low-frequency dynamics. A At that is too large will reduce the 
method’s ability to detect small variations in pilot model parameters. In 
the current study, the length of the time window is chosen to be 20 s. 

In Reference 4, to increase the likelihood of finding a global optimum 
solution, the initial parameter set for the Gauss-Newton optimization is 
determined by a genetic optimization procedure. This procedure requires 
a considerable amount of computational power, especially when performed 
at every time step. To reduce the computational effort in the current 
study, the initial parameter set was constructed from the simulated pilot 
model parameters at every time step. 



t , S 


Figure 4: MLE time window. 


V. Simulation Setup 

A simulation was set up in Matlab® and Simulink® to allow for an evaluation of the estimation of 
time-varying pilot model parameters using wavelets and MLE. This section describes the control task, time- 
varying dynamics, and forcing function and remnant signals used in the simulation. Values for time-varying 
parameters used in the simulation are not directly taken from previous experiment data. However, the 
parameters were adjusted manually to approximate pilot performance results of previous studies. 

A. Control Task 

Figure 5 depicts the simulated closed-loop target-following control task. The task of the pilot is to minimize 
the error e - that is, the difference between the target forcing function f t and the system output y - visible 
on a compensatory display. The time-varying vehicle and pilot dynamics are given by H c (s,t) and H p (s,t), 
respectively. The pilot control signal u serves as the input to the controlled vehicle dynamics and is the 
summation of the linear pilot response function output and a remnant signal n that accounts for non-linear 
behavior. 

The time-varying vehicle dynamics are given by: 
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pilot 



Figure 5: Simulated closed- loop control task with time- varying dynamics. 


H c (s,t) 


K c {t ) 

T c (t)s 2 + s 


(15) 


where K c (t) is a time-varying gain and T c (t) is a time-varying time constant. These dynamics are a single 
integrator below the break frequency of 1/T c (t) rad/s and a double integrator above this frequency. 

The time-varying pilot dynamics are determined by: 


H p {s,t) 


K v (t)(l + T l (t)s)e~ ST ' 


^ nm 

T 2CnmW nmS + s 2 


(16) 


The equalization dynamics of the pilot consist of a time- varying gain K v ( t ) and a time- varying lead time 
constant T/(t). The time delay t v . and the neuromuscular parameters (nm and ui nrn , represent the pilot’s 
limitations and are constant in the current simulation. 

The pilot remnant signal n was simulated by passing a zero-mean Gaussian white noise signal through a 
time-varying low-pass filter: 


H n (s,t) 


K n (t ) 

T n s + 1 


with a time- varying gain K n (t) and a fixed time constant T n . 


(17) 


B. Time- Varying Dynamics 

A single simulation run with a length of 90 s is divided into three parts. In the first part, until 1 1 = 30 s, the 
controlled vehicle and pilot dynamics are constant. Next, until G = 70 s, the dynamics transition from one 
state to the other by linearly changing the time-varying parameters K c (t), T c (t), K v (t), TJ(f), and K n {t). 
After the transition phase the dynamics remain constant again for the remainder of the simulation run. 

The values for the time- varying parameters for the length of a run are given in Figure 6. The fixed values 
for the pilot model time delay, and neuromuscular damping and frequency, are set to t v = 0.20 s, ( nm = 0.20, 
and oj nm = 10 rad/s, respectively. Note that the values for the lead time constants of the vehicle and pilot 
model dynamics are equal throughout the run. This means that a perfect pilot compensation is assumed for 
the double integrator vehicle dynamics at higher frequencies; that is, the open-loop dynamics ( H a i = H p H c ) 
have single integrator characteristics for the entire frequency range. 

Frequency responses of the time-varying vehicle, pilot and open-loop dynamics for every 5.0 s in the 
transition phase are given in Figure 7. As can be observed from Figure 7a, the vehicle dynamics change from 
mostly integrator dynamics to mostly double integrator dynamics. As a result of this transition, the pilot 
needs to generate more lead, as can be observed in Figure 7b. Note that despite the fact that the parameters 
change linearly, the magnitude and phase of the vehicle and pilot responses do not change linearly. Both 
magnitude and phase change more rapidly at the beginning of the transition phase. 


C. Target Forcing Function 

The target forcing function is a sum of sines constructed using the following equation: 
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IH C , deg \H a \, 



Figure 6: Values of the time- varying simulation 
parameters. 


Table 1: Forcing function properties. 


k, - 

n t , - 

Wt, rad s 1 

A t , deg 

(fh, rad 

1 

5 

0.384 

0.629 

-2.571 

2 

13 

0.997 

0.629 

-1.059 

3 

27 

2.071 

0.629 

1.736 

4 

41 

3.145 

0.629 

2.060 

5 

53 

4.065 

0.629 

-2.790 

6 

73 

5.599 

0.063 

-1.221 

7 

103 

7.900 

0.063 

2.020 

8 

139 

10.661 

0.063 

0.127 

9 

194 

14.880 

0.063 

1.483 

10 

229 

17.564 

0.063 

-0.537 


(a) vehicle magnitude 




10' 1 10 ° 10 1 
uj, rad s — 1 


(b) pilot magnitude 


(c) open-loop magnitude 




(e) pilot phase (f) open-loop phase 




Figure 7: Frequency responses of the time-varying vehicle, pilot, and open-loop dynamics from t\ = 30 s 
until <2 = 70 s for every 5.0 s. 
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(a) error signal (b) control signal (c) state signal 





Figure 8: Time traces of the simulation signals. 


N t 

ft(t) = 51 A t( k ) sin [u t (k)t + M k )] ( 18 ) 

k=% 

with Nt = 10 the number of sine waves, and uit, A t and <pt the frequency, amplitude and phase shift of 
the fc th sine wave, respectively. The measurement time used to construct the forcing function was set to 
T m = 81.92 s. With a simulation sampling frequency of 100 Hz, this measurement time contains the highest 
power-of-two data points. The sinusoid frequencies w t (fc) were all integer multiples of the measurement-time 
base frequency, = 2tt /T m = 0.0767 rad/s, and were covering the frequency range of human control (0.1-20 
rad/s). The amplitudes for the lowest five frequencies were set equal and the amplitudes at the remaining 
frequencies were set to one tenth of the low frequency amplitudes. The forcing-function phase distribution 
was randomly generated. The final forcing- function signal was scaled to have a variance of 1.0 deg. Table 1 
provides a summary of all the forcing-function properties. 


D. Pilot Remnant 

The value of the constant remnant filter time constant T n was set to 0.2 s. The time- varying gain K n (t) was 
calculated to induce a certain pre-set power ratio between remnant and pilot control signal (P n = cr^/a^) 
for the entire run. The gain can be calculated using the following expression for the variance of a signal x: 


a x = ~ 
7T 


1 p+oo 1 r+oo 

- / S xx {u)duj = — — / X*(u)X(u})duj 
7T Jo nN Jo 


(19) 


Using this expression to calculate the variances of n and u, and substituting the frequency responses of 
these signals as defined by Equations 20 and 21, the value of K n (t,) can be calculated as a function of the 
time- varying pilot and controlled dynamics. 


N(ju}) = H n (juj)W(juj) 


( 20 ) 


U(ju) 


H P {ju) 

1 + H p (joo)H c (juj) 


F t (ju) + 


Hn(ju) 

1 + H p (ju])H c (ju) 


W(juj) 


( 21 ) 


Time traces of the different signals in the control loop for a single simulation run are given in Figure 8. 
In this simulation run, the remnant was scaled to be 10% of the total pilot control signal (P n = 0.1), using 
the above equations. 


VI. Results 

This section provides the results calculated using the simulation discussed in Sec. V. Simulations with and 
without pilot remnant signals were performed to investigate the influence of pilot remnant on the estimates 
of pilot control behavior and the pilot model parameters. 
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(a) WVL pilot magnitude 


(b) MLE pilot magnitude 





Figure 9: Time-varying pilot identification results from the wavelet and MLE methods (P n = 0). 
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A. Pilot Control Behavior 

Figure 9 depicts the time- varying pilot frequency responses identified using the continuous wavelet transform 
and MLE. For the wavelet calculations, only the last 81.92 s of a simulation run were used; that is, the largest 
power of two data points. The MLE procedure was performed every two seconds, starting at t = 12.0 s until 
t. = 80.0 s. Note that a time trace of 20.00 s was used for the estimation of the pilot frequency response at 
each time step. The time-domain data used for the identification was taken from a simulation where no pilot 
remnant was simulated (P n = 0). The figure clearly indicates the time- varying nature of the pilot frequency 
response. 

For the frequency response identified using wavelets, some edge effects can be observed at t = 0.0 s and 
t = 81.92 s in both the magnitude and phase plots, Figures 9a and 9c, respectively. In addition, some 
oscillatory behavior of the wavelet transform can be observed in time at higher frequencies. Note that the 
time-varying wavelet frequency response is only defined at the input frequencies of the forcing function. 

The magnitude and phase of the pilot frequency responses identified using MLE are given in Figures 9b 
and 9d, respectively. The frequency response appears to be smoother, as knowledge of the pilot model is 
incorporated into the estimate. By comparing the wavelet and MLE magnitude responses around t = 30.0 
s, it can be observed that the rapid change in the pilot frequency response at the start of the parameter 
transition is not captured very accurately using MLE. The maximum likelihood estimate shows a more 
gradual change due to the sliding time window. The sliding time window effectively produces an average 
frequency response over the length of the time window. At the end of the transition phase at t = 70.0 s, 
where the rate of change in dynamics is more gradual, the MLE procedure is capable of accurately estimating 
the change in dynamics. 

To better assess the accuracy of the wavelet and MLE estimates and to allow for a better comparison 
between the two methods, the estimation results are given at three time instances in Figure 10. The time 
instances coincide with the start, middle, and end points of the parameter transition phase. The pilot 
frequency response data is equivalent to the data presented in Figure 9. It can be observed that the pilot 
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(a) pilot magnitude at t = 30 s 


(b) pilot magnitude at t = 50 s 


(c) pilot magnitude at t = 70 s 





(d) pilot phase at t = 30 s (e) pilot phase at t — 50 s 




(f) pilot phase at t = 70 s 



Figure 10: Pilot identification results at different time instances (P n = 0). 


frequency response estimates from the wavelet transform - only defined at the input frequencies of the forcing 
function - follow the simulated pilot response very closely at the given time instances. The estimate from the 
MLE procedure at t = 30.0 s does not follow the simulated response very accurately. This is the same effect 
as observed in Figure 9. The MLE frequency response estimates at the other two time instances accurately 
follow the simulated responses. Note that the data is from a simulation without simulated remnant. 

Figure 11 depicts the estimated frequency responses for a simulation with added pilot remnant. The 
variance of the remnant signal n was scaled to 5% of the pilot control signal variance (P n =0.05). The 
remnant percentage was kept constant over the length of the simulation run using the procedure described 
in Sec. V.D. When comparing the wavelet and MLE ferquency responses, it can be observed that the wavelet 
estimates are most heavily affected by the addition of a remnant signal. A remnant variance percentage of 
5% of the pilot control signal is relatively low. When averaging data runs from an experiment with human 
subjects, the average remnant variance percentage is above 10%. 

B. Crossover Frequencies and Phase Margins 

Figure 12 provides the estimated time-varying crossover frequencies and phase margins from the open-loop 
responses calculated with the pilot frequency response estimates from Figure 9. The values calculated from 
the simulated pilot and controlled dynamics are given for reference. For the calculation of the crossover 
frequency - where the magnitude of the open loop is 1.0 - a single interpolation is required. For the 
determination of the phase margin - the phase difference with 180 deg at the crossover frequency - another 
interpolation is required. The errors between the estimated and the simulated crossover frequencies and phase 
margins are a result of inaccuracies from the frequency response estimates and the interpolation procedures. 
Note that the MLE frequency response estimates have a much higher resolution in the frequency domain 
compared to the wavelet estimates, reducing the interpolation errors. 

The crossover frequency and phase margin estimates resulting from the wavelet frequency responses 
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(a) WVL pilot magnitude 


(b) MLE pilot magnitude 





Figure 11: Time- varying pilot identification results from the wavelet and MLE methods (P n =0.05). 
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closely follow the simulated values. The edge effects of the wavelet estimates in Figure 9 also slightly affect 
the wavelet estimates of the crossover frequency and phase margin at the edges. The crossover frequencies 
and phase margins resulting from the MLE data show a large discrepancy around t = 30.0 s. This coincides 
with the error in the estimated frequency responses observed in Figures 9 and 10. At t = 70.0 s where the 
change in dynamics is more gradual, the values are estimated correctly. 

The estimated crossover frequencies and phase margins from the simulation data with pilot remnant are 
given in Figure 13. This figure again shows that the data resulting from the wavelet estimates is mostly 
affected by the addition of pilot remnant. The error in the MLE frequency response estimates at t = 30.0 s 
is still present. 

C. Parameter Estimates 

The parameters of the pilot model structure defined in Eq. (16) were estimated using the wavelet pilot 
frequency response estimates and the cost function defined in Eq. (9). The parameters were estimated at 
a two second interval to reduce the computational effort. Note that this is a second step to be performed 
after the calculation of the frequency responses using wavelets. The MLE procedure directly estimates the 
pilot model parameters from the time-domain data. For both parameter estimation techniques, the initial 
parameter vector at every time instance was chosen to be the simulated set of parameters. The wavelet and 
MLE parameter estimation results for all five parameters of the pilot model are given in Figure 14. These 
results are for a simulation without simulated pilot remnant. 

From Figure 14 it can be observed that both techniques produce parameter estimation results that 
follow the general trends of the simulated parameters. The pilot gain resulting from the wavelet estimation 
procedure is estimated too low for almost the entire run. The MLE parameter estimates generally fluctuate 
around the simulated values. At the end of the simulation run, around t = 80.0 s, the MLE estimates of 
pilot gain and lead time constant start to deviate from the real value. 

The parameter estimation results for a run with simulated remnant is given in Figure 15. The remnant 
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(a) crossover frequency 



(b) phase margin 
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Figure 12: Estimated time-varying open-loop characteristics from wavelet and MLE methods (P„= 0). 


(a) crossover frequency (b) phase margin 




Figure 13: Estimated time- varying open-loop characteristics from wavelet and MLE methods (P„=0.05). 


variance percentage was kept constant at 5% for the length of the simulation run. It can be observed that the 
estimates from both the wavelet and MLE procedure are affected by the remnant. The wavelet parameter 
estimates are mostly affected by the pilot remnant, especially before the start of the parameter transition 
(t < 30 s). The lead time constant estimated using the two-step wavelet procedure starts to deviate from 
the middle of the run. The pilot gain and lead time constant resulting from the MLE procedure start to 
deviate from the end of the transition phase. 

Simulation runs were performed with percentages of pilot remnant variance up to 20%. However, no 
parameter estimates could be obtained using the two-step wavelet parameter estimation procedure. The 
bias and variance of the pilot frequency responses calculated using the wavelet transform proved to be too 
high to produce reliable parameter estimates in the second step. In most cases the constrained optimization 
procedure returned parameter values at the boundaries of the parameter search space. The MLE procedure 
was able to produce reliable parameter estimation results comparable to the results given in Figure 15. 

VII. Discussion 

Simulations of a closed-loop manual control task with time- varying dynamics were performed to evaluate 
a two-step parameter estimation technique using wavelets and a one-step parameter estimation technique 
using MLE for the estimation of the time-varying pilot model parameters. The simulation was set up using 
knowledge of pilot performance and equalization characteristics. However, no parameter values from real 
experiments were taken. The use of simulations allowed for an objective comparison of both parameter 
estimation procedures and their sensitivity to pilot remnant. Future studies should investigate the use of 
both techniques in a human-in-the-loop experiment. 
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(b) lead time constant 
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Figure 14: Estimated time- varying pilot model parameters from wavelet and MLE methods (P„=0). 


The pilot remnant was simulated as a low-pass filtered Gaussian white noise signal. The time constant 
of the filter was kept constant. The variance percentage of the remnant in the pilot control signal was 
also constant during a simulation run. These values were kept constant to reduce the amount of variables 
changing at the same time to allow for a better evaluation of performance in estimating the time-varying 
pilot equalization parameters. However, in a real-life task, the characteristics of pilot remnant would change 
over time in correspondence to the changing vehicle dynamics. For example, the percentage of pilot remnant 
is generally higher when controlling double integrator dynamics (acceleration control) compared to single 
integrator dynamics (velocity control). 

In the case of no simulated pilot remnant, the wavelet parameter estimation technique was able to provide 
accurate estimates of the time- varying pilot frequency response, and crossover frequencies and phase margins. 
The results were only slightly affected by the edge effects that are commonly seen when using wavelets. 8 
However, the frequency response functions identified using wavelets were greatly affected with the addition 
of simulated pilot remnant, even with a variance as low as 5% of the total pilot control signal variance. 
The parameter estimates resulting from the two-step wavelet parameter estimation technique proved to be 
unreliable in this case. With higher variances of pilot remnant, no parameter estimates could be obtained. 

The performance of the wavelet parameter estimation method is dependent on the type of wavelet function 
and the selected wavelet parameters. The Morlet wavelet used in this study is the most commonly used 
wavelet function for the identification of human control behavior. 9, 10 However, depending on the task, other 
types of wavelet functions may provide better results. In addition, time averaging the wavelet transforms 
before the calculation of the pilot frequency response functions might improve the method’s accuracy when 
applied to data with considerable pilot remnant. 8 The averaging of the wavelet transform was not performed 
in the current study. 

Without the addition of pilot remnant to the simulated control task, the MLE procedure provided 
accurate results for the pilot frequency responses and model parameters. However, the method proved to be 
less effective when the dynamics of the pilot are changing rapidly, as seen by the significant error between 
simulated and estimated response functions at t = 30.0 s in Figure 10. This is caused by the averaging effect 
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Figure 15: Estimated time-varying pilot model parameters from wavelet and 


MLE methods (P n =0.05). 


of the moving time window. In general, the pilot frequency response and model parameter estimates from 
the MLE procedure were less affected by pilot remnant. This warrants the use of the MLE procedure for 
future studies concerning the estimation of time-varying pilot model parameters. However, more knowledge 
needs to be gained on the effects of the moving time-window length on the estimation results. 

In the current study, only an unconstrained gradient-based Gauss-Newton optimization was performed 
to find the parameter estimate with the highest likelihood using MLE. In Reference 4, it was found that a 
genetic algorithm in addition to the Gauss-Newton optimization - significantly improved the MLE method’s 
capability of finding a global optimum solution. However, the genetic algorithm significantly increases the 
computational effort needed, especially when performed every time step. A solution could be to use the 
genetic algorithm to calculate an initial parameter vector at the start of a run only and use the results of the 
Gauss-Newton optimization as the initial parameter vector for the Gauss-Newton optimization in the next 
step. However, if the optimization at a single time instance would produce a non-optimal solution, it might 
not be possible to achieve global optimum solutions at future time steps. 

To improve the performance of MLE for the identification of time-varying parameters in future studies, 
the method should be enhanced to estimate the parameters as a function of time directly as opposed to 
using a sliding time window. This would reduce the computational effort and increase accuracy. However, 
assumptions need to be made on the parameter functions and properties before the analysis. These assump- 
tions can easily be made in an experiment where the controlled time-varying dynamics are known, but this 
will be more difficult in real-life applications. 

Using both estimation techniques, the accuracy of the estimation results are highly dependent on the 
values of the pilot model parameters in combination with the frequency content of the forcing function. 
For example, when the inverse of the lead time constant of the pilot model decreases to a value below the 
lowest input frequency of the forcing function, the pilot model gain can be estimated less accurately. This 
problem can be minimized by careful design of the forcing functions and prior knowledge of the change in 
time-varying dynamics. 

The current study was a first attempt to evaluate wavelets and MLE for the estimation of time-varying 
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pilot model parameters. More research into the performance of both methods is needed. In addition, the 
methods need to be extended to a multi-loop case; that is, pilot models with different perceptual inputs. 
In addition, the variation of parameters other than the pilot equalization parameters, such as the pilot 
time delay and neuromuscular parameters, should be considered in future studies. Finally, the variation in 
parameters was assumed to be linear in the current study. Nonlinear variations in parameters should be 
considered in future research. 

The estimation of a time-varying component of pilot control behavior parameters adds a very complex 
problem to the already difficult problem of pilot control behavior identification. Future experiments on 
time-varying pilot control behavior should be designed carefully. The accuracy of the methods discussed 
here should be evaluated more thoroughly before using the methods in experiments or real-life task analysis. 

VIII. Conclusions 

This paper presents two methods for the estimation of time-varying pilot model parameters, a two- 
step method using wavelets and a windowed maximum likelihood estimation method. The methods were 
evaluated using simulations of a closed-loop control task with time-varying pilot equalization and vehicle 
dynamics. Simulations were performed with and without remnant. Both methods give accurate results 
when no pilot remnant is present. The wavelet transform is very sensitive to measurement noise, resulting 
in inaccurate parameter estimates when considerable remnant is present. The MLE method is less sensitive 
to pilot remnant, but is unable to detect fast changes in pilot control behavior. 
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